Gauge cooling in complex Langevin for QCD with heavy quarks 
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We employ a new method, "gauge cooling", to stabilize complex Langevin simulations of QCD with 
heavy quarks. The results are checked against results obtained with reweigthing; we find agreement 
within the estimated errors. The method allows us to go to previously unaccessible high densities. 
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Introduction. - An important problem of particle 
physics is to derive the phase diagram of hot and dense 
QCD from first principles. The difficulty lies in the fact 
that the action becomes complex at finite density. Vari- 
ous methods have been employed to get around this prob- 
lem: Taylor expansions at /i = and analytic continua- 
tion from imaginary /i, reweighting etc., but these have 
a limited range of applicability A general method 
that recently has received a lot of attention is the Com- 
plex Langevin Equation (CLE) 0-13 . It has had some 
successes as well as some failures; the failures seem to be 
related to the problem that the system wants to spend 
"too much time too far out" in the complexified configu- 
ration space; this problem and possible ways to deal with 
it were discussed in @, |^ . 

In QCD at finite density the enlarged gauge freedom in 
the complexified field space can actually open a new way 
to limit the dangerous large excursions by employing non- 
unitary gauge transformations. This freedom was used in 
a somehwat different way for a SU{2) toy model [9|; the 
procedure used there, however, does not generalize to a 
lattice model in an obvious way. 

Here we introduce the method of gauge cooling (g.c), 
designed to keep the link variables as close as possible to 
the unitaries, thereby preventing the dangerous large ex- 
cursions. We test the idea first for a simple Polyakov loop 
model, where we have exact results to compare our data 
with; then we apply it to a heavy quark approximation 
of QCD with chemical potential (HQCD) (cf. 0). 

Complex Langevin for gauge models - The complex 
Langevin method for a complex action S is based 

on setting up a stochastic process on the complexification 
of the configuration space. The longtime average of holo- 
morphic observables is then supposed to give the correct 
average corresponding to the complex weight exp(— 5). 
In lattice models of QCD the configuration space is 
5.[;(3)#ii„ks (ggg for the Real Langevin approach); 
after complexification this becomes iS'L(3, C)'^''"'''^ [13| . 



The complex Langevin equation (CLE) is 

idU^,f,)U-l = - ^ XaiDa,x,f.Sdt + dWa,x,^) , (1) 
a 

Aq, a — 1 . . . , 8 are the Gell-Mann matrices; dwx,^,a are 
independent Wiener increments, normalized as 

{dWa,x,ti{t)dWa',x',ti'{t')) ^ '^Saa' Sw Sxx' S f,^,' dt ; (2) 

Da,x,n is a differential operator (derivation) acting on 
functions / of S'C/(3)#""'"'^ as 

Da,x^,f{{U}) = hm - [f{{U{6)}) f{{U})] , (3) 

where {U{5y\ means the variable Ux,^ has been replaced 
by ex.'p{i5Xa)Ux,fi with all other variables unchanged. 
The first term on the rhs of Eq.([T]) (the drift term) is 
gauge covariant and transverse to the gauge orbits since 
S is gauge invariant, but the noise term contains compo- 
nents along the gauge orbits. Since the gauge group is 
noncompact, the process (1) may go far from the unitary 
submanifold. 

An Euler discretization of Eq.([T]) (see also [6j) gives 

Ux,i^ exp |-^iAa(eii'a,a:,,i + Ve^ya.x./i)! Ux,^j. , (4) 

where Ka^x,^, = Da,x,iiS is the drift force and rj are inde- 
pendent Gaussian noises satisfying 



{"n^jX ,x' '^^aa'^xx'^^ 



(5) 
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Gauge cooling. - It is a well known problem of the CLE 
method that the system may drift too far out into the 
imaginary directions, causing numerical problems. To 
quantify the distance from unitarity we use the 'unitarity 
norm' 

F{{U}) ^Y.^r[Ul^Ux,, + {Ul.r'U-}, - 2] > , 

(6) 

which vanishes if and only if all the f/'s are unitary. 
Gauge invariance is the freedom to transform variables 

Ux,f. ^ v-^Ux,^yx+[. (7) 
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for all links x,x + fl, where initially are uni- 

tary, but after complexifieation t^, V^+/i £ SL{3,C). We 
use this freedom to minimize the unitarity norm, with- 
out changing the observables. Concretely, we intersperse 
standard 'dynamical' Langevin sweeps with several 'g.c. 
sweeps'. These are deterministic moves essentially in the 
direction opposite to the gradient of F. 

A maximally nonunitary gauge transformation at site 
y in direction a is given by 



Uy-fi,i_t n- Uy-f,^f_, exp(-Q!Aa) 



(8) 



(a g R) for all fj,, with all other link variables remaining 
unchanged. The 'gauge gradient' of F in the a direction 
at the lattice site y is then 



Ga.y = Da^yF — 2trAa 

+ 2trAa 



Uy.fiUy^^ Uy_p^^^Uy-iX^fi 



The g.c. updates of the configuration are given by 



(9) 



Ux.~L i-^- exp 



Ux- 



Ux-ji,,i exp ^ aXaG 



(10) 



where a — ea; a determines the strength of the g.c. force, 
whereas e is a discretization parameter as in Eq. (j4|) . Note 
that even if a is not small, Ea. ([TU| is still a gauge trans- 
formation; it just might not be optimal for reducing F . 

Polyakov loop model. - The Polyakov loop model is 
given by a ID lattice consisting of N links with periodic 
boundary conditions. Analytically it reduces to a one- 
link integral, but it is a useful laboratory to check the 
effect of g.c. The action is given by 



5 = /3itrC/i...[/jv, +/32trC/^' 



,.c/r 



(11) 



where we allow (3i^2 to be complex. Here we choose /3i = 
/? -I- Ke'^, 1^2 = (3* + Ke^^; S will in general be complex. 
The observables of interest are trT^*-' = tr(J7i . . . JJjv)'^, 
k — ±1, ±2, ±3. The effect of the g.c. on (trT') is shown 
in Fig[l] where we vary a of Eg. dTOl) . In Fig. [2] we show 
the effect of varying a on the distribution of the values 
of trV. We see that with too small a the distribution 
has a 'skirt' of slow decay. From [8] we know that slow 
decay typically leads to incorrect results; a small skirt, 
however, does not lead to appreciable deviations from 
the exact values. The main point is that enough cooling 
leads to correct results, see FiglT] 

We also measured tr [/^^ and tr U^^ with similar re- 
sults, for N up to 1024; it turns out that increasing N 
requires an increase in cooling. 
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FIG. 1: Polyakov chain; average Polyakov loop vs. a. 
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FIG. 2: Polyakov chain: histograms of Polyakov loops. 

Heavy quark QCD (HQCD). - The model was first ex- 
plored using the CLE without g.c. in [7] on small lattices 
and for moderate /i. It is defined by dropping all spatial 
hopping terms of the quarks; this amounts to simplify- 
ing the Wilson fermion determinant of lattice QCD with 
nonzero chemical potential /i to 

detM(/x) EE \Y Det(I -f- C7'^)^Det(l + C'V'^f , 

X 

(12) 

where C = [2Kexp(^))]^', C = [2k exp(-^)]^' , Det 
refers to the color degrees of freedom and 



Vx 



r=0 



HQCD is described by the action 

5= §5g ({£/})+ In det M(/i) 
6 



(13) 



(14) 



where Sg = J2p^^ tr(L/p) is the Wilson plaquette ac- 
tion and antiperiodic b.c. in time are chosen for the 
fermions. We have 

Det (I + CVxf = (1 + + iCPx + iC^P'x? ■ (15) 
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FIG. 3: HQCD: evolution of the unitarity norm with and FIG. 4; HQCD: Baryon density and average phase for /3 



without g.c. for /3 — 5.9, 8x6 lattice. 



5.9, 8'' X 6 lattice. 



and similarly for the second factor in Eg. dT^ . where 



P. 



(16) 



A related model has been studied numerically by a 
reweighting (RW) technique to deal with the complex 
density in (and more recently in [IBl), where for 
larger fj, and larger lattices the sign problem causes bad 
signal to noise ratios. See [13] for the general setting. 

We simulate the model by the CLE on various lattices, 
with K = 0.12, various values of /3 and fi, using Eq.Q 
with e = 10~^ — 2 X 10"''. In some cases we use adap- 
tive control for the dynamical stepsize and the number 
of cooling steps. The system is thermalized after a cold 
start up to Langevin time t — 10 before taking averages. 

In Figl3]we show the Langevin evolution of the unitar- 
ity norm: it is seen how g.c. stabilizes the process. We 
see that g.c. is needed even for /i = where the process 
is real, but has an instability pushing it away from the 
unitary submanifold in the absence of cooling. 

In Fig U] we show the baryon density n/usat as a func- 
tion of the chemical potential /i on an 8^ x 6 lattice. We 
find saturation for /i > 2, suggesting that our results are 
correct also for fairly large ^. The figure also shows the 
average phase of exp(— S*), defined by 



dctM(^) 



(17) 



^dctM(-/i) , 

(see 0])- With growing density the average phase drops 
to almost zero and rises again in the saturation region. 

FigEl shows (P) and (P') vs. /i. Both reach a max- 
imum where the density takes off (/i ~ 1-4), then drop 
again to nearly zero (see also [16l - [l8| ): the curve for {P') 
is shifted slightly to the left with respect to (P). The 
analytic curves for /3 = in the figure exhibit a similar 
behavior; they are calculated essentially as in [l^. 

This behavior is due to the fact that for small as well as 
large fj, the determinants in Eqs. (|12p are dominated by 
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FIG. 5: HQCD: (V), {V) vs. at ^ = 5.9 on a 8^ x 6 lattice; 
solid lines: analytic strong coupling result. 



the constant terms, leading to small expectation values of 
P and P'. In physics terms, to get a nonzero answer one 
has to form a localized colorless bound state between the 
dynamical quarks and the external charge provided by P 
or P'. For small fi > 0, when typically there is only 1 
quark present, it can form a 'meson' with P' but not with 
P; for larger fi, when typically 2 quarks are present, P, 
but not P' can form a 'baryon' with them. The situation 
is slightly complicated by the fact that up to 6 quarks 
can exist at every site, so P' can also form a colorless 
state with 4 and P with 5 quarks. With no quarks or 
in the completely filled state it is not possible to form 
such bound states. The simple picture given here ignores 
fluctuations, including a small probability, damped by ^, 
of finding antiquarks, so at /i = (P) and (P') are small 
but nonzero, whereas at very large fi they will go to zero. 
Thus both quantities will have a peak, but (P) will peak 
later than (P'). 

Comparison of CLE with RW. - To check the CLE 
results we also simulated HQCD with a RW technique, 
in the region where this is feasible (cf. [l3|)- We com- 
pare the results for the Polyakov loops and for plaquette 
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FIG. 6: HQCD: Comparison of (P), (P') and plaquettes for 
RW and CLE at /3 = 5.9, 6* lattice, a = 1, 12 g.c. steps. 
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FIG. 7: HQCD: (P), (P') and plaquettes vs. P, ^^ = 0.85 on 
a 6"^ lattice with adaptive step size. 



expectations in Figl6] The results agree within the er- 
rors (large errors only occur for RW). Note that around 



fi = 1.2 the RW method breaks down, whereas CLE with 
g.c. works fine for arbitrarily large fi. 

In FiglT] we show the expectation values of Polyakov 
loops at fi — 0.85 vs. /3, both obtained from RW and 
CLE. The figure shows that CLE works for large (3 and 
also allows us to cross over into the confining region. 
Deeper in the confining region, however, we see small 
differences between the RW and CLE data, apparently 
related to an emergence of a 'skirt' of the distribution. 

Concluding remarks. - The instabilities that arise 
without or with too weak g.c. are apparently caused by 
(1) the existence of repulsive fixed points, pushing the 
configurations exponentially fast away from the unitary 
submanifold and (2) roundoff and other numerical errors. 
Each gauge orbit contains configurations arbitrarily far 
from the unitaries, and there roundoff errors pile up to 
affect also the observables, violating gauge invariance. 

The modified process remedies these problems; it 
thereby also avoids the slow decay of the equilibrium dis- 
tribution. The CLE method with g.c. is tested success- 
fully for the Polyakov loop model. It allows, apparently 
for the first time, to simulate a real QCD like gauge model 
at finite density all the way into the saturation region. 
The method does not suffer from any sign or overlap 
problem. There is no fundamental obstacle to extend 
the method to fuU QCD. 
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